% 读取CSV文件，路径中的反斜杠需要转义
filename = 'C:\\Users\\王鹏翔\\Desktop\\实验4\\math-practice\\4\\temp.csv';
data = dlmread(filename);

% 获取数据维度
[n_days, n_hours] = size(data);

% 构造设计矩阵X和响应向量T_vector
X = [];
T_vector = [];
for d = 0:n_days-1       % 天数从0到25（共26天）
    for h = 0:n_hours-1  % 小时从0到23
        % 计算周期性项
        angle = 2 * pi * h / 24;
        cos_term = cos(angle);
        sin_term = sin(angle);

        % 构建设计矩阵行
        X_row = [d, 1, cos_term, sin_term];
        X = [X; X_row];

        % 获取对应温度值（注意Matlab索引从1开始）
        T_val = data(d+1, h+1);
        T_vector = [T_vector; T_val];
    end
end

% 线性回归求解参数
theta = X \ T_vector;
k = theta(1);       % 季节趋势斜率
c = theta(2);       % 合并的截距项 (b + a0)
a1 = theta(3);      % 余弦项系数
a2 = theta(4);      % 正弦项系数

% 计算拟合值
T_fit = X * theta;

% 绘制实际温度与拟合温度对比图
t = 1:length(T_vector);
figure;
plot(t, T_vector, 'b.', t, T_fit, 'r-', 'LineWidth', 1.5);
xlabel('时间点（每小时数据）');
ylabel('温度 (°C)');
title('沈阳气温实际值与模型拟合对比');
legend('实际温度', '拟合温度', 'Location', 'Best');
grid on;

% 绘制季节性趋势和周期性波动分解
seasonal_trend = X(:,1)*k + c;       % 季节性线性部分
periodic_effect = T_fit - seasonal_trend; % 周期性波动部分

figure;
subplot(2,1,1);
plot(t, seasonal_trend, 'g-', 'LineWidth', 1.5);
title('季节性趋势 (线性部分)');
ylabel('温度 (°C)');
grid on;

subplot(2,1,2);
plot(t, periodic_effect, 'm-', 'LineWidth', 1.5);
title('日周期性波动 (三角函数部分)');
xlabel('时间点（每小时数据）');
ylabel('温度 (°C)');
grid on;

% 输出拟合参数
disp('拟合参数:');
disp(['k (季节趋势斜率) = ', num2str(k)]);
disp(['c (截距 b+a0) = ', num2str(c)]);
disp(['a1 (余弦系数) = ', num2str(a1)]);
disp(['a2 (正弦系数) = ', num2str(a2)]);
